function idat = fsb_normalize_baseline_add(idat,sandbox,hrf_pred)

% FSB - DEV :  Additive baseline normalization
%
% EXAMPLE:
% idat = fsb_normalize_baseline_add(idat, sandbox,2)
%
% INPUT:
% idat:     4D-image information
% sandbox:  sandbox struct with experiment information
%
% OUTPUT:
% idat:    normalized image information
%
% CALLED BY:
% FSB.m
%
% METHOD:
% calculates the mean baseline intensity of every voxel in every trial by
% taking all the datapoints at which the first regressor was still less than
% a 20th of the maximum response in every trial
% calculates the mean baseline of all trials
% calculates the difference between the mean baseline for every trial and
% the mean baseline for all trials. The result is a vector for every voxel
% and every trial
% Adds this difference to the voxel intensity values, thereby achieving a
% normalization of mean trial values.
%
% NOTES:
% By commenting the first loop, now the trial baselines are compared to the
% mean brain intensity instead of the mean baseline intensity. Thus, the
% scaling does essentially scale up the whole intensity, as in an ideal
% case the baseline should usually be lower than the later response.
%
% Copyright 2010 MPI for Biological Cybernetics
% Author: Steffen Stoewer
% License:GNU GPL, no express or implied warranties
% 
% $Revision 1.0
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

h = waitbar(0,'Normalizing trials with own baseline intensity...');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Transform int16 to single for further processing
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
idat = single(idat);
idat(idat==0)=NaN;
idat(idat<=0)=NaN;
blength=zeros(1,max(sandbox.intrial(:,3)));
vxcorr = zeros(1,max(sandbox.intrial(:,3)));

%~~~~~~~~~~~~~~~~
% Calculate mean baseline and trial value at every voxel
%~~~~~~~~~~~~~~~~
hrf_max = max(sandbox.hemodynamics(:,hrf_pred));
bsl_ind1 = find(sandbox.hemodynamics(:,hrf_pred)<hrf_max/20);
bsl_ind2 = find(sandbox.intrial(:,2)<max(sandbox.intrial(:,2))-3);
bsl_ind = intersect(bsl_ind1,bsl_ind2);
bsl_mean = nanmean(idat(:,:,:,bsl_ind),4);

for x = 1:max(sandbox.intrial(:,3));

    waitbar(x/max(sandbox.intrial(:,3)));
    scanind = find (sandbox.intrial(:,3)==x);
    
    if min(scanind)==0;
        scanind(1)=[];
    end

    idat_red = idat(:,:,:,scanind); % Gives volumes in trials
    bsl_trial = sum(sandbox.hemodynamics(scanind,1),2); % Doing baseline calculation based on first regressor time course

    % Looks for baseline volumes by looking for values that are below
    % 1/20th of maximum response
    bsl_max = max(bsl_trial);
    bsl_ind = find(bsl_trial<(bsl_max/20));
    bsl_ind(bsl_ind>6) = [];
    
    bsl_ind2 = mean(idat_red(:,:,:,bsl_ind),4);

    voxel_corr = bsl_mean-bsl_ind2; % correction factor voxel by voxel
    voxel_corr(voxel_corr==inf)=1;

    voxel_corr_4D = repmat(voxel_corr,[1 1 1 size(idat_red,4)]);
    idat_red = idat_red+ voxel_corr_4D;
    idat(:,:,:,scanind) = idat_red;

    vxcorr(x) = nanmean(voxel_corr(:));
    blength(x)= size(bsl_ind,1);

end;

idat = int16(idat);
disp(['Mean baseline length: ' num2str(mean(blength))]);
disp(['mean voxel corr: ' num2str(nanmean(vxcorr))]);
close(h);

vxcorr = vxcorr/nanmean(idat(:))*100;
figure; bar(vxcorr);
set(gcf,'Name','Trial baseline intensity correction');
xlabel('Trial');
ylabel('Average signal deviation %');

disp(' ')
end
